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Five numerical methods were compared for calculating two-dimensional, transient natural con- 
vection in an enclosure. Both implicit and explicit procedures were considered. Requirements for 
numerical stability were derived from analysis and experience, and when satisfied, the calculated flows 
for all methods were found to be similar. Consideration was also given to the accuracy and (energy and 
vorticity) conservation of the methods. One method was found to be conservative and stable without a 
restriction on the spatial mesh increment. This method can be successfully applied to nonlinear flows, 
but care must be exercised due to the presence of truncation errors which introduce false transport 
mechanisms. 
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1. Introduction 
1.1. General 

The purpose of this paper is to compare the practical stability, truncation errors, and conser- 
vation properties of five finite-difference procedures for solving a problem in natural convection. 
Three of these methods are from the current literature [2, 6, 20J l and two are developed here. The 
ability of the various methods to produce physically meaningful solutions is examined; and solutions 
calculated with the various methods are compared. 

The physical problem chosen for study is the natural convection flow induced in a vertical 
circular cylinder by a small hot spot centrally located on the floor (see fig. 1). Transient and steady 
laminar flows are considered in two-dimensional, axisymmetric cylindrical coordinates. Such a 
fluid motion is described by three simultaneous partial differential equations [4, 10]: (a) an equation 
relating stream function to vorticity, eq (2); (b) a time-dependent equation for vorticity, eq (3); 
and (c) a time-dependent equation for temperature, eq (4). The first equation is of elliptic type, and 
the last two are of parabolic type. 

1.2. Numerical Methods Tested 

This section provides a brief introduction to the five numerical methods tested. All methods 
employ successive over-relaxation [16, see chapter 11] for solving the finite-difference approxima- 
tion of the elliptic equation, and no difficulties are encountered. The parabolic equations, on the 
other hand, contain nonlinear terms with first-order derivatives which express the influence of 
convection (hereafter called the convection terms). These terms introduce serious problems of 
stability and conservation into the finite-difference scheme. The numerical methods differ only in 
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their treatment of the parabolic equations, and these differences are highlighted in table 1. Methods 
IV and V are new and are modifications of methods II and I, respectively. 

TABLE 1. Numerical methods applied to parabolic equations for temperature and vorticity 

Checks mean that the numerical methods possess the indicated properties. 



Method 


Reference 


Explicit 

or 
implicit 


Difference form for convection terms 


Order 
of trun- 
cation 
error 


No time 
step re- 
striction 


No spatial 
mesh size 
restriction 


No false 
diffusion 
and con- 
vection 


Satisfies 
conser- 
vation 


I 




Barakat and Clark [2] 


E 

1 
E 

1 
E 


2 pt forward or backward 


h 

h 2 
h 2 
h 2 
h 




\S 






II 


Wilkes and Churchill [20J 

Fromm [6] 




l* 
V 
i* 




III 








^ 


IV 


Present study 








^ 


V 


Present study 


3 pt noncentral forward or backward.. 




" 


\^ 













An inspection of table 1 reveals that no one method is conservative, free of false transport, 
and stable without a restriction on the spatial mesh size. For conservation reasons, methods III— V 
are to be preferred over methods I and II. Of these, methods III and IV require a reduction in 
spatial mesh size with increasing flow velocities to achieve stability, while method V does not. 
This mesh size restriction leads to prohibitively large demands on computer time and storage. 
The stability of method V, together with its conservation, are essential for calculating flows where 
the convection terms are important, as in the present problem. The introduction of false diffusion 
and convection, however, means that the calculated flows must be interpreted with care. In general, 
this false transport is important only in those flow regions where methods III and IV were found 
to break down. 

Since the completion of this work, two additional differencing schemes have come to the 
author's attention [15, 1]. The first of these [15] was independently developed for steady flows and 
employs a differencing of the convection terms identical to method V. The second method [1] has 
been employed for meteorological problems and was developed for the vorticity equation. This 
method conserves vorticity, mean kinetic energy, and mean square vorticity. The convection terms 
are differenced to h 2 accuracy, as with methods II- IV in table 1, but conservation of mean square 
vorticity should lead to improved stability characteristics. 

A comparison of the methods listed in table 1 is presented in the following sections. The physi- 
cal problem is posed in section 2, some features and properties of the numerical methods are 
discussed in section 3, and comparative physical results are presented in section 4. 

2. Mathematical Description of the Physical Problem 

Consider the motion of a viscous fluid within a vertical circular cylinder of height a and radius 
b (see fig. 1). Erect a cylindrical coordinate system (x, r) with origin at the center of the base. The 
flow is assumed to be axisymmetric with no variations in the azimuthal direction. The fluid is initially 
motionless and at a uniform temperature TV The enclosure walls are also at this temperature, 
except for a small centrally located circular spot on the base of radius c which is at a temperature 
T h > T . The temperature difference T h -T {) initiates and sustains the natural convection flow within 
the enclosure. 
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FIGURE 1. Cylindrical enclosure and coordinate system. 
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The Boussinesq approximation is used [4, see p. 16); in this density (p) is assumed constant 
except for the generation of buoyancy forces. Other fluid properties are taken as constant: kine- 
matic viscosity (*>), thermal diffusivity (k), and volume expansion coefficient {(3). We introduce the 
following dimensionless quantities: time, r=(/cA/ 2 )/; vertical and radial coordinates, A ' = x\a and 
R = rja\ vertical and radial components of velocity, U={gIk)u and V=((iIk)v; and temperature 
e=(T-To)l(T h -To). 



The governing equations in dimensionless form are 

l a^ i a^ 

~R aft' R dX' 



i a 2 ^ a (i dv\ 

ft dX 2 dR \R aft/' 



aft , [d(Uil) , a(Fft) 



dr 



-+ 



n 



)X 



+ ■ 



dR 



GrPrt^+Pr 



dm a /i d(Rsi) 

dX 2 dR \R dR 



(1) 

(2) 
(3) 



and 



ae 



+ 



a(i/0) l a (ft KB) 

dX R aft 



a 2 e l a /„ ae 
a* 2 ft aft \ aft 



(4) 



Equation (3) contains the Prandtl number Pr=vJK and the Grashof number Gr= z g^(Th~~ To) a 3 /p 2 . 
The latter number is used for natural convection flows and denotes the product of buoyancy and 
inertia forces divided by the square of the viscous force. The acceleration of gravity is denoted 
by g. The existence of a stream function ^V is assumed such that the velocities are given by (1). 
The mass conservation equation div U=0 is then automatically satisfied. The vorticity vector 
11= curl U has only an azimuthal component, ft = BV/dX— dU/dR, given by (2). The conservation 
equations for vorticity and energy are (3) and (4) respectively. In these equations, bracketed terms 
on the left and right sides respectively denote the convective and diffusive transport terms. The 
term containing Gr in (3) represents the vorticity source due to buoyancy. 

The convective terms in eqs (3) and (4) are in a "conservation" form [3, 6] suitable for methods 
III, IV, and V of table 1. Alternate "nonconservation" forms for these terms are utilized for methods 
I and II. The alternate equations are obtained by introducing a modified vorticity ft' = (l/ft)ft into 
(2) and (3), followed by subtraction of ft' div U=0 and 9 div U= from (3) and (4) respectively. 
In terms of modified vorticity ft', eqs (1) through (4) become 



1 a^ 
j_ 2JL V-- 

ft aft ' 

l a 2 ^ 



l a¥ 

R dX' 

m \r c>r 



a') 

(2') 



an' n dCi' { v dQ,'_ 

Bt dX dR 



dQ 



R dR 



+ Pr 



dm' i d 

dX 2 R dR 



G 



i d(Rm') 



dR 



(3') 



and 



dr dX dR " dX 2 



L±( R §9\ 

R dR \ dR/' 



(4') 



Equations (l)-(4) and (l')-(4') are subject to initial conditions and boundary conditions. The 
initial conditions are: 



a = 9 = 0for t<0, OssA'ssl, andO=sK ^ R b . 
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(5), (5') 



The boundary conditions 2 for t 3= are: 



<fr=d<iridX=o 

e=i 

9 = 0.5 

e=o 
v=dvidX=e=o 



for* = 0, 

for* = 0, 
forA-0, 
forJT = 0, 
forZ=l, 

forfi = 0, 
for R = ft ft , 



« R =s /?„; 
0s/{< R c ; 

R = R r ; 

tf, < ft =s R„: 

all/?; 

all X; and 

all*. 



(6), (6') 



For equations (l')-(4') replace fl= in (5) by 12'= 0, and 0= (R = 0, all A') in (6) by the dual con- 
ditions V=dU/dR = 0. a 

The boundary conditions introduce two geometric parameters: the aspect ratio of the enclosure 
(radius/height) R b = bin and the relative size of the heat source (heat source radius/enclosure height) 
R c — c/a. 

To determine and ft' along the centerline, special forms of the conservation equations are 
needed to avoid an indeterminate form as R — » 0. By incorporating boundary conditions at R = 
and using L'Hospital's rule, eqs (4), (3'), and (4') respectively reduce to 



ae awe) 

3t dX 



3 HRVO) 
~ dR* 



a 2 9 d^e 
dX* dR*' 



^r +u w— GrPr w +Pr 



a 2 ft' em' 

ax* DR* 



and 



(It dX 



d 2 G a 2 e 

dX* BR*' 



(7) 

(T a) 

(7'b) 



For convenience, the velocities U and V are retained explicitly in the problem formulation. 
However, an examination of eqs (l)-(7) (or alternately, eqs (l')-(7')) reveals that the essential 
dependent variables are ^, A, and (or ^P, A', and 0), while the independent variables are 
X, R, and r. The parameters of the problem are R&, R c , Pr, and Gr. Throughout this study, the 
aspect ratio, relative heat source size, and Prandtl number are held fixed at R&=1, /? c = 0.1 and 
Pr = 0.7, respectively. The Grashof number Gr is assumed equal to 1 X 10 5 , except in section 4.4, 
which considers larger values of this parameter. 

3. Numerical Methods, Formulation and Some Properties 
3.1. Grid System and Calculation Sequence 

An approximate solution of eqs (1)— (7) or (l')-(7') will be obtained at a finite number of 
grid points having coordinates X=iAX, R=jAR, and at discrete times r n , where i 9 j, and n are 
integers. The grid spacings in the X and R directions are denoted by hX and A/?. The symbol 
T n denotes the time level after the nth time step Ar w . The values of *P, fl, 0, U, and V at each 
grid point should be thought of as average values over a small volume of fluid surrounding the 
point. 



2 Note that explicit boundary conditions for vorticity on the solid boundaries are not available. This causes no difficulty with the explicit solution methods and only 
slight difficulty with the implicit methods. 

3 The centerline boundary condition £l = K£Y = is satisfied for all finite values of fl'. Although the vorticity fl is zero, the modified vorticity fl' is not, in general, 
zero. Barakat and Clark [2] used the modified vorticity ft' and assumed ft' — along the centerline. This is incorrect, and may have contributed to the oscillations in 
stream function (near the centerline) and heat transfer observed in their calculations. The vorticity ft' on the centerline can be evaluated with the aid of eq (7' a). 
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All of the numerical methods listed in table 1 advance the fields of ^P, (1, and across a time 
step in the same sequence. All quantities are assumed known at a time r n . The sequence for ad- 
vancement from time T n to the new level t w+ i = t w + At is as follows: 

1. Temperature (B) at all interior grid points is advanced with suitable difference approxima- 
tions of (4) and (7) or (4') and (7'b). This is discussed in section 3.2. 

2. Vorticity (fl or ft') at all interior grid points is similarly advanced with approximations of 
(3) or (3') and (7' a). This is discussed in section 3.2. 

3. Stream function (fy) at all interior grid points is brought up to date with the new vorticity 
field by using a difference approximation of (2) or (2'). This is discussed in section 3.3. 

4. The vorticity on the solid boundaries is determined from the new stream function field, as 
are the velocities U and V. This is discussed in section 3.3. 

5. The field of mesh points is scanned in order to determine the time step At for the next time 
advancement. This is limited by practical stability requirements and is discussed in section 3.4. 

3.2. Parabolic Equations 
a. Differencing Schemes 

The approximations of the equations for temperature and vorticity (for steps 1 and 2 above) 
are described in conventional notation in table 2. In this section we comment on the distinguishing 
features and some of the properties of these schemes. 



Table 2. Finite-difference approximations of the parabolic equations for temperature and vorticity 



ae ao a*e 



<K) d(ue) a*o 



Illustrated with the simplified temperature equations, h U — — — — (for methods 1 and II) and 1 — — (for methods 1 1 1 — V ) . 

dr dX dX* ()t HX dX 2 



ao 



ae a(fyo) 



a 2 Q 
aA 2 



e.fv-e," 



— + 0(At) 



,,J AY ,,J f \u\d*e 



an. ,-2ef,+e? , 



At/2 



■^ + 0(At) 



V — -c7» i - iLLj L -^-f-0(AV) 2 

dX ' ,J 2AA 



B," . ,-20" . + 6" , , 

_JJ L_U + 0(ZLX)2 



(AA) 2 



At/2 



ao o; 1 ;' ,-0?v, 

U — =U? i - Lth2 L ^ + 0(AA) 2 + 

dX ' J 2AA \AV/ 



Vaa/ 



o? + v J -20|' i y + e|'_v.j 



(AA) 2 



+ 0(AA) 2 -r()(^- 



Ill 



er,y-er; 



£ + 0(At) 2 



b(uq) Uf +u pf +iti -uu,Pt-uj 



Bf +UJ - e?y - efo 1 + Of. 



+ 0(AA) 2 



(AA) 2 



+ <aa' )2 + °(at) 2 - 



Same as method 11a. 



d(UB) _Uf +t j9f +UJ -Uf_ ui Qf_ ui 

8X 2AA 



Same as method 11a. 



Same as method Ilh. 



a ( ue ) t/« + , .o» + v , - */»_ , o,»v j / at \ 

IF- - M "-- + o ( ^ + o (-) 



Same as method Ilh. 



Same as method I. 



a((70) 
dX 



AA 
AA 



! , (Ofj, Uf_ UJ <0) 

+ 0(AA) + 0(AY) 2 



Same as method I 



*</„,„ -(</ mil .« + </ m . n )/2. 



Application of the various numerical methods to the parabolic equations is illustrated with 
the one-space dimension temperature equation 



285 



^ + ~lx ajp~° (8) 

in place of eqs (3), (4), and (7); and 

^ +f/ l¥-aF =0 (8) 

in place of eqs (3'), (4'), and (7'). Equation (8) is employed for methods III-V and eq (8') for methods 
I and II (see table 2). It is understood that the discussion and differencing of each term of these 
equations applies to similar terms in the two-space dimension temperature and vorticity equations. 

Methods I, III, and V are explicit procedures for calculating quantities at time n+1 in terms 
of their values at times n and n — I. Methods II and IV are two step implicit procedures (see [12, 
16, 20]), going from time n to n . + 1/2 in step a, and from n + 1/2 to ra-fl in step b. During the 
first half time step (or second half time step), the X and R space derivatives are respectively ap- 
proximated at time levels n and n+ 1/2 (or n+1 and n + 1/2). Clearly, the time steps a and b are 
implicit in the R and X directions, respectively. This is called an alternating direction, implicit 
scheme and leads to a tridiagonal matrix of unknown temperatures or vorticities for each /?-row 
or ^-column, as the case may be. These matrices are readily inverted by a simple algorithm. 

Column two of table 2 lists the approximation of the derivative dO/dr. A central time differ- 
ence is used with method III, and forward time differences with all others. The last column of the 
table lists the approximation of the diffusion terms, d 2 0/dX 2 . Three point space differences are 
used for all methods. 

The approximations of the convection terms UdQ/dX or d(UO)/dX are listed in column three. 
At any grid point the velocities are evaluated at time n and are treated as constants over the time 
step. With method I, the convection terms are approximated with two-point forward or backward 
differences as the coefficient velocity U is positive or negative, respectively. Methods II— IV all 
employ three-point central differences for these terms. Method V employs a modified form of 
forward or backward differences in which the mean velocity U n , m in the numerator is multiplied 
by On, m or B/,+1, m as U n , m is positive or negative, respectively. If the two mean velocity coefficients 
are of different sign, one term from the numerator of each of the two approximations shown is 
required. 

Two points about the vorticity equation should be noted. First, this equation contains a buoy- 
ancy source term which is not listed in table 2. The term contains dQ/dR, which is approximated 
with three-point central differences. Method III evaluates the term at time n, all other methods 
evaluate it at time /i + 1. Second, when the implicit methods II and IV are applied to the vorticity 
equation, it is necessary to temporarily assume the vorticity distribution on the solid walls at 
times n+1 and n + 1/2 equal to that at time n. Thus, the wall vorticity is out of step with the ad- 
vancement of the interior vorticity field. The error introduced by this assumption decreases with 
both small time steps and the approach to steady state. 

The rest of this subsection presents brief comments on the stability and conservation prop- 
erties of the parabolic difference equations. Additional discussion is provided in sections 3.4 and 
4.3, respectively. 

The stability of the various numerical methods is principally due to some key approximations 
shown in table 2. For methods I and V, it is the use of forward or backward differencing for the 
convection terms [2, 7, 14 (see p. 194)]. For method III, it is the duFort-Frankel differencing [5] 
of the diffusion terms d 2 Q/dX 2 . For methods II and IV, it is the alternating direction, implicit 
nature of the scheme (due to Peaceman-Rachford [12, 16 (see p. 366), 20]). 

A lengthy but straightforward study of methods I-V shows that conservation is satisfied for 
III-V, but not for I and II. The elements of such a study are presented in [3, 6]. Conservation of 
energy or vorticity within the grid system exists if the difference equations for temperature or 
vorticity are summed over all interior grid points and no spurious sources or sinks of these quan- 
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tities are found. The convection and diffusion terms cancel in pairs. The net energy or vorticity 
transport from the wall mesh points into the enclosure just balances the net increase of energy 
or vorticity within the mesh system. Heat transfer rates to and from the enclosure can thus be 
determined by calculating the transfer from wall mesh points. If conservation is satisfied, there 
is no net increase of energy within the enclosure at steady state, and the total heat added to the 
enclosure balances the heat removed. 



b. Truncation Errors 

The truncation errors shown in table 2 are worthy of some discussion. These estimates are 
obtained by Taylor series analysis as described in references [5, 11, 13]. Methods II— IV are most 
easily discussed and are treated in the first paragraph. Methods I and V are more complicated, 
and are discussed in the second and third paragraphs, respectively. 

For methods II— IV, all space derivatives are approximated to (AX) 2 accuracy. The use of 
multiple time levels in conjunction with these space derivatives, however, introduces additional 
truncation errors of 0{Ar/AX) for methods II and IV and of 0(Ar/AX) 2 for method III. Substitution 
of the difference approximations listed in table 2 into eqs (8) or (8'), as appropriate, introduces 
truncation errors of 0(Ar)+ 0( At/A*) + (MM) 2 for methods II and IV and of 0(Ar) 2 + 0(Ar/A#) 2 
H-O(AA') 2 for method III. The coefficients of these terms involve derivatives of higher order than 
those of the basic eqs (8) or (8'). Stability requires that At <^ AJ(, and when AA' is small, these trunca- 
tion errors are usually neglected. 

With method I, the O(AA') approximation of the convection term UdO/dX introduces a large 
truncation error. A Taylor series expansion for 0f-ij about 0,,j can be rearranged to solve for 

(deidXhj 



(-) ■ 

\dXjij 



9;,j — ©;- 



i.j 



AX 



+ 



/^B\ (AX)* (d 3 G\ 1 



ax /d-e 

2 



(9) 



The first term on the right side is a backward difference. The term in braces is the truncation 
error. Multiplication of eq (9) by Uij(Uij> 0) leads to the form shown in the table. The solution 
of the difference equation (i.e., the approximation of (8')) is equivalent to a solution of the differ- 
ential equation 



dO f[ dO_ 

bt ax 



ll/U 



U\\d 2 G 



2 7 l)X 2 



:()(At)+0(AA'K 



(10) 



The coefficients of the terms 0(At) and Q{AX) 2 involve derivatives of higher order than appear 
in (8'), and for reasons noted in the preceding paragraph, these truncation errors are usually 
neglected. The differencing of the convection terms in equation (8') introduces an additional, 

Jyl Only for small AX ^ 



or false, heat diffusivity AA' 1 ~. Only for small AX -—-* does the difference approximation ap- 



proach the differential equation (8'). For large AX — , care must be exercised because a false 

heat diffusivity is introduced (or a false viscosity in the case of the vorticity equation). 4 

Turning now to method V, the analysis is similar to that presented for method I above. The 
truncation error of O(AA^) in the convection term d(UG)/dX is more complicated, however. It can 
be shown that the solution of the difference equation (i.e., the approximation of (8)) is equivalent to 
a solution of the differential equation 



dr 



+ 



[ a me ) ax dUdO) 
1 ax ~ 2 axax]' 



+ A* 



Ul± 



dX*~ 



0(At) + 0(A*) 2 



(11) 



•This problem cannot bt- overcome by eliminating d-O/dA' 2 from (9) with, for example, three-point central differences. If this is done.eq (9) becomes a three-point 
central difference approximation of dO/dX, with its associated stability problems. 
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when Ujj, £/j_i, j > 0. Again, the terms on the right side are usually neglected. There is now a 

false heat diffusion (or a false viscosity) as well as a false convection. The latter is a subtraction 

AX 
from the real convection and is small whenever -r- (dU/dX)(dOldX) is small. When this last term 

and AX — are small, the difference approximation approaches the differential equation (8). When 
these terms are large, care must be exercised. 

3.3. Elliptic Equation, Velocities, and Wall Vorticity 

The computations described in steps 3 and 4 of section 3.1 are identical for all five numeri- 
cal methods, and are integrated and presented in this section. Equations for the vorticity fl (used 
for methods III— V) are considered here. Appropriate equations in terms of modified vorticity 
II' (for methods I and II) are obtained by the substitution fl = Ril' . Step 3 is considered first. 

The elliptic equation (2) relates stream function to vorticity at any given time. The new vortici- 
ties ftf/5 1 are introduced and the space derivatives are approximated by three-point central differ- 
erences with truncation errors of O(AJf) 2 and Q(AR) 2 . An iterative technique known as successive 
over-relaxation [21] was employed to obtain the new stream function field. No problems of stability or 
conservation are encountered. Thus, if ^^ denotes the approximate stream function at a point 
interior to the solid boundaries and centerline after s iterations, a further approximation ^ r ^j 1) 
is obtained from 



*<->= a -^fS+ _ L _ _ j _V z _ j\ {**«& 1+ m 

(Axr + (ARr\j+t + j-i) 



x i+i,j x f-i,j 



+ l 



(A/0 2 






7^2 J 2 



(12) 



It is understood that all values of ^P in this equation pertain to time n+l. The optimum value of 
the relaxation parameter co was calculated for a given system of grid points (see eqs (1.9) and (4.6) 
of reference [21]). 5 

Successive sweeps of eq (12) over the field of mesh points were alternated in the vertical and 
radial directions. Iteration was terminated when 

max i^-^lli 

-^ < 0.0001 . (13) 

max i^:; 1 '! 

hj 

The number 0.0001 was arrived at by experimentation, noting that more stringent criteria led to 
no essential difference in the solution of test problems. The numerator on the left side of (13) is 
the maximum change in stream function occurring in the field of grid points as a result of one ap- 
proximation step. The denominator is the maximum value of stream function in the field. 

In practice, the most time consuming part of the whole calculation was the iteration of eq (12) 
at each time step. One iterative sweep took almost as long as the entire computation of the new 
temperature or vorticity field described in the previous section. The following procedures resulted 
in a considerable savings in computer time. A maximum number of iterations, s max , was permitted. 
For 100, 400, and 1000 mesh points, s max was 20, 25, and 35, respectively. The convergence criterion 
(13) was generally achieved within s max iterations after the initial flow transient was over. The 
stream function field was then updated less frequently, typically every 1+0.25 (s max -s conve rgence) 
time steps. 



5 An independent determination of the optimum a> was made with the aid of a computer. Equation (12) was solved for various to starting with a given ft field. The 
value of o) leading to fastest convergence was sufficiently close to the calculated value so that the latter was used. 

288 



The computations described in step 4 (sec. 3.1) are discussed in the rest of this section. The 
numerical approximations up to this point, together with the boundary conditions (6), provide the 
entire fields of ^ and O, and the field of O interior to the solid boundaries, all at time n + 1. Still 
to be calculated at time n + 1 are the unknown values of vorticity II on the solid boundaries 6 and 
the velocity fields U and V. The new stream function field is used for this purpose. 

The vorticity on the solid boundaries is obtained by first reducing eq (2) to Q = — (llR)d 2J ¥/dX 2 
for X = or X= 1, and to fi== — (l/Rb)d 2j ¥ldR 2 for R = Rt>. The second order derivatives are approxi- 
mated by expanding ^ in a Taylor series about the wall. Expansions for the two grid points nearest 
the wall, together with the boundary conditions on ty, yield wall vorticity approximations such as the 
following, which applies whenX= 0: 

Finally, the fields of U and V are calculated with three point central difference approximations 
of eqs (1). These approximations (of 0(A/Y) 2 and O(Aft) 2 ) are not presented here, but it is readily 
verified that they automatically satisfy a central difference approximation of the continuity equation, 
div = 0. The conservation of mass within the grid system is thereby established. 

3.4. Practical Stability Considerations for the Parabolic Equations 

At this point, the fields of ^P, ft, G, £/, and V are current at time n+l for all methods. Further 
integration of the parabolic temperature and vorticity equations in time and space requires con- 
sideration of practical stability. This is step 5 of section 3.1. Practical stability imposes a restriction 
on the size of the time step At for methods I and V (eqs (19) and (20), respectively) and restrictions 
on Ar, A/V, and A/? for methods II- IV (eqs (15)). The size of At is calculated from the appropriate 
equations (assuming for methods II- IV that the restrictions on AX and AR are satisfied). The fields 
of ty, ft, B, £/, and V are then advanced across this time step by repeating the whole cycle described 
in sections 3.2 and 3.3. 

A complete analysis of the nonlinear equations to determine the exact form of the stability 
requirements is not always possible. Such is the case for methods II— IV. The linearized stability 
analysis of von Neumann [11] has been applied to methods II [19] and III [6] and leads only to a 
restriction on At for method III. Numerical experience with methods II— IV for a limited range of 
flows suggests that the following empirical restrictions be applied at each grid point for integration 
of the temperature and vorticity equations: 

M ~TO' M ~ra' (15a) 



At* 
and 



_2_ + _2_ 

|_(A¥) 2 (Aft 



(15b) 



lAY 1 



1/7- -I \V- 
AX AR ' { } 



Note the restriction on spatial mesh size (15a), which presumably results from the use of three 
point central differences for the three methods. Within limited ranges these restrictions have led 
to stability, but they cannot be regarded as general. In particular, during flow transients with large 
initial changes, time steps smaller than suggested by (15b) or (15c) were required. 



8 This delay in calculating the wall vorticities causes no difficulty with the explicit methods I, III, and V. For the implicit methods II and IV, it was necessary in 
section 3.2a to temporarily assume that the wall vorticities at times n + 1 and n+ 1/2 were equal to their values at time n. 
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Analytical stability criteria are possible for methods I and V. The analysis is applicable to 
explicit, nonlinear, two time level differencing schemes and is discussed in reference [2]. This 
paragraph is based on that analysis. To be specific, stability of the integration of the temperature 
eqs (4) or (4) is considered here; that for vorticity (3) or (3) is similar. In methods I and V, the un- 
known temperature 0f ? + 1 can be written as an explicit linear combination of computed values at 
time n: 



e?y = ai6f +1}i + a*e?_ x tj + a 3 e? fj + a 4 e/; i+1 + ^0? 



i-i 



(16) 



where the a* denote coefficients which vary in time, but which are constant over a time step. 
The essential feature of the method is to require that the norm of the matrix of coefficients au at 
all times is bounded by unity. The row norm is used, and leads to 



max \ 2 



dk\ 



1 



(17) 



where max denotes the maximum value of the quantity in braces for all i,j in the grid system. Sta- 
i,j 

bility in the sense of Lax and Richtmyer [9, 14 (see p. 44)] exists if inequality (17) is satisfied. For 
methods I and V, satisfying (17) is equivalent to requiring that all the coefficients a* are positive. 
This follows because for methods I and V the coefficient a^ has the special form: 



fl.j = 1 — {a i ■+• a 2 + a± -+- «5 } 



(18) 



which can be readily verified. The conditions under which the coefficients au are positive for 
methods I and V are discussed next. 

For method I, the coefficients (ik (k= 1, 2, 4, 5) are always positive, whereas a% can be made 
positive by restricting the size of the time step At. For Pr ^ 1, the most severe restrictions on At 
are: 

; + -ttt^I , (19a) 
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and 



At = 
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AX (AXf (ARf\ 



(19c) 



These restrictions respectively follow by requiring that a.j be positive for the finite-difference 
representations of eqs (4'), (7'b), and (7' a). 

For method V (as for method I above), the coefficients cik (A=l, 2, 4, 5) are always positive. 
In order to make as positive, however, several different restrictions on At arise because method V 
employs forms for the convection terms that depend upon the sign of the mean velocities, as noted 
in table 2. For the case when the mean velocities £/;,j, £Ai,j, Vjj and Vij-i are positive 7 and 
Pr ^ 1, the greatest restrictions on At result from the difference forms of the temperature equations 
(4) and (7): 



A^ 






+ 1 + 



1\ V x 
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7 The mean velocities are defined by 



Um, n = and Vm, n ~ ~ 
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and 



At-- 



^0+4^ + -^ + ^ 
AX AR (AX)* (AR) 2 ] 



(20b) 



Additional forms similar to (20) arise when the mean velocities are negative, or a combination of 
positive and negative. 

Some observations on the application of the stability requirements can be noted. It is imprac- 
tical to alter AX and AR according to eqs (15a) during a computation to achieve stability for methods 
II- IV. Accordingly, a cut and try procedure is employed. The time step, however, is adjusted during 
a computation with all methods, using equations (15b), (15c), (19), or (20), as appropriate. The field 
of mesh points is scanned in order to determine the largest allowable At which will satisfy the 
stability requirements. Time steps smaller than this value can be used, of course, to provide a 
safety factor. In practice, 80 percent of the allowable At was used with methods I-IV and 95 percent 
with method V. 

4. Comparative Physical Results 

Natural convection flows calculated with the five numerical methods are presented and com- 
pared in sections 4.1 through 4.3 for a Crashof number of Gr= 1 X \(Y\ The calculations employed 
a uniform grid spacing of AX = AR = 0.05, with a total of 121 mesh points. The calculation of physi- 
cally meaningful flows at higher values of Gr, and some questions of convergence, are considered 
in section 4.4. 

Flow patterns and temperature fields are illustrated in figures 2 through 6 with computer- 
drawn graphs displaying sets of streamlines and isotherms. The location of these streamlines and 
isotherms was determined by linear interpolation of the computed values at the mesh points. In 
each of the graphs, the centerline of the cylindrical enclosure is shown on the left. The abscissa 
is the radial coordinate /?, and the ordinate is the axial coordinate X. The heat source on the floor 
is denoted with a thick line between R = and ft = 0.1. 

4.1 . Transient Streamline and Temperature Fields 

All the numerical calculations were carried from the initial quiescent condition forward in 
time until steady state was achieved. The transient and steady-state flows calculated during this 
process are qualitatively similar for all five numerical methods. The purpose of this section is to 
discuss the physical nature of the flows, consequently, only the transient results for one method 
need be considered. Method V is considered representative, and the transient flow and temper- 
ature fields calculated with this method are illustrated in figures 2 and 3, respectively. These figures 
are each composed of four graphs, arranged in order of increasing time, t. 
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FIGURE 2. Transient streamline fields calculated with method V at various times (r), Gr— I X 10 5 . 

Streamline evolution is qualitatively similar for all five numerical methods. The walls and centerline correspond to ^ = 0; the dot has value W miix ; the remaining 
streamlines correspond successively to ^ values of 0.1, 0.3, 0.5, 0.7, and 0.9, of Vmw 
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FIGURE 3. Transient temperature fields calculated with method V at various times (r), Gr= 1 X 10 5 . 

Temperature evolution is qualitatively similar for all five numerical methods. Temperature (G) is the curve parameter; the heat souree (thick line on floor) is 
at - 1 . 

Immediately after the start of heating, a ring vortex of warm fluid forms near the origin (figs. 
2a and 3a). This vortex rises and moves radially outward from the centerline (figs. 2b and 2c). A 
gradual development into the steady state flow and temperature fields then follows. The steady 
state fields reveal a ring vortex centered just above ^ = 0.5, R = 0.5 (fig. 2d) which is driven by the 
heated fluid which rises along the centerline from the floor to the ceiling (fig. 3d). 

A sensitive indication of the approach to steady state was provided by the overall energy bal- 
ance on the enclosure. When the heat transfer rates into and out of the enclosure were within 1 
percent of an asymptotic value, steady conditions prevailed. Steady state was achieved by r=0.25 
with all methods. 

4.2. Comparison of Steady-State Streamline and Temperature Fields 

The qualitative similarity of flows calculated with the five methods was mentioned in the last 
section. The purpose of the present section is to illustrate this similarity in a quantitative way, by 
comparing steady state flows obtained with methods I-V. Such a comparison of the steady-state 
streamline and temperature fields is presented in figures 4 and 5, respectively. The sets of stream- 
lines and isotherms shown have the same numerical values as in figures 2 and 3. Figures 4 and 5 
each contain four individual graphs, which respectively pertain to numerical methods I, II, III, 
and V. Computed results from methods III and IV agreed to at least three significant figures, and 
graphs of the steady state streamline and temperature fields for the two methods are identical. 
Consequently, only the results of method III are shown. Figures 4d and 5d duplicate figures 2d 
and 3d respectively. 
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FIGURE 4. Steady-state streamline fields for the various numerical methods, Gr= 1 X 10 5 . 

Results for methods III and IV are identieal. The walls and centerline correspond to ty= 0; the dot has value ^ max : the remaining streamlines correspond suc- 
cessively to ^ values of 0.1, 0.3, 0.5, 0.7, and 0.9, of ^ max . 
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FIGURE 5. Steady state temperature fields for the various numerical methods, Gr= 1 X 10 5 . 

Results for methods III and IV are identical. Temperature (O) is the curve parameter; the heat source (thick line on floor) is at = 1. 

A comparison of figures 4a-4d or 5a-5d reveals a striking similarity in the streamline and 
temperature fields computed with the various numerical methods. Method I shows the greatest 
departure from the group with the largest value of ^ max (caption, fig. 4a) and in the shape of the 
isotherms above the heat source (fig. 5a). The physical reason for this is that method I transfers 
more heat by convection from the heat source than the other methods (and thus has the greatest 
rate of fluid circulation, ^P^x)- This is discussed further at the end of section 4.3. 

Figures 4 and 5 illustrate flows in which both convective and diffusive transport exists. The 
contribution of convection to the temperature fields in figure 5 can be appraised by referring to 
figure 6. This figure illustrates the steady state temperature field for static conduction (no fluid 
motion) and was obtained with method V by setting Gr= U-- V— ft = ^ = in the basic differential 
equations (l)-(4). The conduction fields calculated in this way for numerical methods I-V are 
identical (to at least four significant figures) because all methods approximate the heat diffusion 
(or conduction) terms to the same spatial accuracy. 




05 X 



FIGURE 6. Static conduction temperature field at steady state, all methods. 

Temperature (Ol is the curve parameter; the heat source (thick line on floor) is at O = 1 . 



Table 3 lists the computer time (in seconds) required to calculate the flows illustrated in figures 
4 and 5. The particular computer used had execution times of 1.75 and 2.625 /xs for nine digit float- 
ing point addition and multiplication, respectively. The only significant observation that can be 
made about this table is that method III required the least computer time, and method I the most. 
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Differences between methods could be reduced by altering the safety factor used for the time step 
At or by more efficient programming. 

TABLE 3. Computer time required for test case with Gr= 1 X 10 5 and a 21 X 21 mesh 



Method 


CoiT 


iputer time 


I , 




sec 
335 


II 

Ill 

IV 


270 
230 
309 


V 


285 







4.3. Heat Transfer Rates 

The heat transfer rates calculated by the five numerical methods are conveniently discussed 
in terms of the rate of heat addition (Q m ) or heat removal (Q out ) from the enclosure. These are total 
quantities, obtained by integration over the areas of the heat source and cold walls, respectively. 
The integration is performed by computing the heat transfer by convection and conduction from 
wall mesh points to adjacent mesh points, as discussed in the paragraph at the end of section 3.2a. 
A dimensionless heating rate can be defined as 

where X is the thermal conductivity of air, a is the height of the enclosure, and AT is the imposed 
temperature difference driving the flow, AT=T h — T {) . Equation (21) defines a dimensionless heat 
addition (3> in ) or heat removal (4> out ) as the rate of heat transfer is () in or Q oui , respectively. 

The heat transfer rates discussed in this section pertain to a ramp temperature change at the 
edge of the heat source. This ramp is the grid approximation of the boundary conditions (6) (a step 
change) and appears as a linear variation of temperature along the floor from = 1 at R — Ri, — AR/2 
to O = at R = Rt>+ AR /2. Both numerically and experimentally, a step change in temperature is 
difficult to achieve. The rate of heat transfer, while being finite for a ramp, would be infinite for a 
step. Twenty-one radial mesh points were used throughout this study and an extension of this 
study [18]. The corresponding ramp (see fig. 4 of [18]) is a close approximation of the floor tempera- 
ture profile in a physical experiment [17] with which the numerical flows could be compared. 

Curves illustrating the heat transfer rates for numerical methods I and II as a function of time 
are presented in figure 7; similar results for methods III— V are presented in figure 8. Figures 7 and 8 
illustrate nonconservation and conservation methods, respectively. The ordinate is the dimension- 
less heat transfer rate into (4> in ) or out of (0 0Ut ) the enclosure. The abscissa is the dimensionless 
time r. 

The calculated heat transfer rates due to static conduction (Gr=Q) are virtually identical for 
all five numerical methods. These conduction curves are shown in both figures 7 and 8 as a basis 
for comparison. Note that the heat removed (dashed line) approaches the heat added (solid line) 
with increasing time. Thus, the finite-difference approximation of the diffusion terms for all methods 
conserves energy within the grid system. 

Methods I and II, figure 7, employ difference approximations of the convection terms which 
do not conserve energy. For method I, the rate of heat removal exceeds the rate of heat addition 
at steady state, indicating a net production of heat within the enclosure. For method II, the heat 
addition exceeds the heat removal, indicating a net absorption of heat within the enclosure. Clearly, 
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Figure 7. Dimensionless heat flow into (4>j„) and out of ($> oll t) the enclosure as a function of time (r) for methods I and 

II (Gr=i X 1(P) and static conduction {all methods). 



methods J and II violate the laws of thermodynamics. These differences in the rates of heat addition 
and heat removal are due to residues which remain when the difference equations for temperature 
are summed over all interior grid points, and the convection terms do not cancel in pairs. Methods 
I and II thus do not retain the energy conservation expressed by the basic differential equations. 8 

Methods III through V, figure 8, employ difference approximations of the convection terms 
which conserve energy. Clearly, for each method the rate of heat addition equals the rate of heat 
removal at steady state. Methods III and IV employ derivative approximations with similar spatial 
accuracy. As a consequence, the heat transfer curves for these methods are virtually identical. 

Two general observations about figures 7 and 8 can be made. First, the shape of the heat 
addition curves is qualitatively similar for all methods, as is the shape of the heat removal curves. 
All curves for heat addition (solid lines) follow the conduction curve initially, then break away and 
achieve a steady value by t= 0.05. The heat removal curves (dashed lines) also follow the conduction 
curve initially. The heat removal rate then falls below the conduction value as it becomes more 
difficult for heat to be conducted radially from the heat source to the nearby cold floor counter to 
the incoming flow. At about r=0.04, the rate of heat removal begins to increase sharply as the 



•The authors of methods I |2| and II |20| did not calculate the heat transport from node to node within the mesh system. Instead, lour point difference formula 

were employed to approximate the temperature gradient at the walls, and thus the rates of heat transfer. For unsymmetric heating, this approach does not lead ti 
equality between the steady state rates of heat removal and heat addition. 
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FIGURE 8. Dimensionless heat flow into (<& in ) and out of (<& (m t) the enclosure as a function of time (t) for methods 111, IV 

and V (Gr = 7 X 10 h ) and static conduction {all methods). 

heated fluid reaches the ceiling above the hot spot. As the fluid spreads radially outward along the 
ceiling, the rate of heat removal continues to rise. The lower curve then gradually approaches the 
upper (in fig. 8). 

The second general observation that can be made about figures 7 and 8 concerns the magnitude 
of the steady-state heat additions. Clearly, the difference between the values for the five methods 
at Gr= 1 X 10 5 and the value for conduction is a measure of the convective heat transfer from the 
heat source. This convective contribution decreases for the five methods in the sequence: method 
I > II > V > 111= IV. The relative magnitudes of these contributions can be estimated by applying 
the difference approximation of the convection terms UdO/dX or d(UQ)/dX (table 2) to the grid 
point (i=l, 7 = 0). A portion of this difference approximation is associated with the convective 
transport of energy from grid point (i = 0, j=Q) to (i = 1, j = 0). A careful study reveals that the 
convective heat transport for the five methods is proportional to the following quantities: method I 
ocf/i,o; II ocf/i, (l + Oi,o)/2; V ccf/ K0 /2; and III and IV oc£/ 1>o lf p/2. (The boundary conditions 
£/o, o = and 0o,o= 1 were employed.) The temperature 0i,o lies between and 1; thus, the afore- 
mentioned sequence of methods is in order of decreasing convection. This agrees with the results 
in figures 7 and 8. 

4.4. Extension to Higher Gr, Convergence 

The calculation of flows at high Grashof numbers is considered in this section. Results at 
Gr = 1 X 10 5 in sections 4.1 and 4.2 revealed a qualitative similarity in the calculated flows for the 
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five numerical methods. Extension to higher Gr, however, means that the methods must continue 
to produce physically meaningful solutions, and this is closely associated with the practical sta- 
bility, truncation errors, and conservation properties discussed in sections 3 and 4.3. For conserva- 
tion reasons alone, we discard methods I and II from further consideration, and consider only the 
conservation methods III— V. 

Of the latter methods, III and IV impose a spatial mesh size restriction to achieve stability 
while V does not. Violation of this restriction leads to oscillations in the fields of 0,^,(7, and V . 
The effect of violating the spatial mesh size restriction is illustrated in figure 9 for the steady-state 
centerline temperature at Gr= 1 X 10 6 . The ordinate is the axial coordinate X and the abscissa is 
the temperature 0. The curves on the right pertain to method III with the number of vertical mesh 
points, M, equal to 21, 41, and 61. (Results for methods III and IV were identical.) For comparison, 
results from method V are shown on the left with M= 21, 51, and 101. In all cases, the number of 
radial mesh points, /V, is 21. 
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FIGURE 9. Steady state distribution of temperature (9) with height (X) along the centerline; methods III and V, Gr= I X 10*. 

Curve parameter is the number of vertical mesh points (M). 



From figure 9, it is immediately apparent that method III leads to an oscillating temperature 
near the ceiling (X= 1) when M = 21. Repeating the calculation with M = 41 reduces the oscillation, 
but a local maximum persists at X = 0.95. Further refinement to M = 6\ eliminates this hot spot. 
With 21 and 41 vertical mesh points, occasional negative temperatures were calculated during 
the transient. The steady state results for M=21 reveal a region of negative temperature atX = 0.8. 
With M=ll, register overflows occurred early in the transient, and a stable integration was not 
possible. 

With method V, no such oscillations or negative temperatures were calculated. Increasing the 
number of vertical mesh points from 21 to 101 reveals that the temperature distribution appears 
to asymptotically approach a limiting curve. Note that the asymptotic temperature distributions 
for methods III and V are not identical; this is attributed to differences in the truncation errors. 

Oscillations such as those in figure 9 are physically unrealistic, and appear only when the mesh 
size restrictions (inequalities (15a)) are not satisfied. Such oscillations have not been observed with 
method V. The inequalities (15a) are satisfied for the test case at Gr= 1 X 10 5 with a 21 X 21 mesh. 
As the Grashof number is increased, it has been found for the present problem [18] that U miix <* 6V 1/2 . 
Clearly, through (15a), this quickly leads to a very fine grid and prohibitively large computer storage 
requirements. 
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The freedom of method V from a spatial mesh size restriction (and from oscillations) sug- 
gests its application to high Gr flows. Due to truncation errors which increase with Gr, such 
results must be interpreted carefully. The solutions should be tested for convergence or compared 
with physical experiment; two such tests for method V are discussed in the subsequent paragraphs. 

A type of convergence check is illustrated in figure 10. The ordinate is the steady state heat 
flow into the enclosure, O in , and the abscissa is the number of vertical mesh points, M. The curve 
parameter is the Grashof number. Solid lines pertain to method V, a single dashed line is shown 
for method III. For reasons noted at the beginning of section 4.3, the number of radial mesh points 
was held fixed at 7V= 21 (except for the one solid data point for which N = '31). The dotted line is an 
analytical result [8] for conduction of heat from a heated disk into a semi-infinite medium. The 
disk has a temperature profile quite close to the ramp profile of the grid approximation. The 
numerical results for static conduction are reasonably close to the analytical conduction value. For 
convective flows at Gr — 1 X 10 5 , 1 X 10 6 , and 4 X 10 7 , the numerical results tend toward an asymp- 
tote with increasing M. 
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Figure 10. Steady state dimensionless heat flow {<Pm) versus number of vertical mesh points (M). 

Number of radial mesh points (N) held fixed at 7V= 21 (except for one solid data point for which /V = 31). 



Additional support for the use of method V is provided by comparison of a numerically cal- 
culated streamline field with physical experiment. An example of such a comparison is presented 
in figure 11 at Gr = 4 X 10 6 . The numerical results are from an extension of this study [18], which 
employed a 51 X 21 grid (M X TV) and covered a range of Gr from 4 X 10 4 to 4 X 10 10 . The experimental 
results are from reference [17] which covered a range of Gr from 8 X 10 5 to 1 X 10 10 . For Gr up to 
about 1 X 10 9 , the comparison between numerical and physical experiment is very good. Above 
this value, the physical experiment indicates turbulence. The numerical experiment was laminar 
(presumably due to false viscosity and the assumed two-dimensional motion) up to Gr = 4X 10 10 . 
At this value, the numerical flow developed a periodic vortex shedding, suggestive of the onset of 
laminar instability. This vortex shedding was adequately resolved with the grid spacing used. 

A significant consideration in the application of any finite-difference procedure to a flow 
problem is the amount of computer time required to reach a steady state. Method V was applied 
in the numerical study of reference [18], and the computing time is shown in figure 12. Clearly, 
as Gr increases, the amount of computing time also increases (as Gr l/4 at large Gr). Approximately 
one-half second was required for each time step, but for large Gr the size of the time step decreases 
rapidly and many more steps are required. 



298 





FIGURE 11. Comparison of numerical results of method I {from reference |18|) with physical experiment {from reference 

[\l])at Cr=4X 10*. 
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FIGURE 12. Computing time required for method V at several Grashof numbers (Gr), numerical study of reference [/#]. 

Calculations carried from quiescent initial conditions to steady state with a V/x/V = 51x21 grid. 
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5. Concluding Remarks 

Five finite difference procedures were compared for calculating transient natural convection 
flows in an enclosure. Although the problem was formulated in axisymmetric cylindrical coordi- 
nates, the conclusions should be applicable to any two-dimensional coordinate system. 

To achieve stability of the time and space integrations, restrictions are imposed on the size 
of the time step for all methods (i.e., I-V) and on the size of the mesh spacing for methods II— IV. 

Energy and vorticity are not conserved within the grid system with methods I and II, but are 
conserved with methods III— V. For this reason, the latter methods are to be preferred. 

Of the conservation methods, method III appears to require less computer time than methods 
IV and V (table 3). Method III is therefore preferable when its associated mesh size restriction can 
be satisfied. If this restriction cannot be satisfied (which is usually the case) oscillations in the 
flow develop. A method free of a mesh size restriction (and oscillations) is then recommended, 
method V. The flows calculated with method V must be interpreted carefully, however, due to the 
truncation errors. 



The author is indebted to L. Orloff, Factory Mutual Research Associate at the NBS, for devel- 
oping a computer program to plot streamlines and isotherms directly from computed numerical 
fields of stream function and temperature. 
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